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We examine in detail the causes of the structural transitions that occur for those small Lennard- 
Jones clusters that have a non-icosahedral global minima. Based on the principles learned from 
these examples we develop a method to construct structural phase diagrams that show in a 
coarse-grained manner how the equilibrium structure of large clusters depends on both size 
and temperature. The method can be augmented to account for anharmonicity and quantum 
effects. Our results illustrate that the vibrational entropy can play a crucial role in determining 
the equilibrium structure of a cluster. 
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I. INTRODUCTION 

The determination of a cluster's equilibrium structure 
is a challenging task for the theoretician. For small clus- 
ters the first aim is usually to determine the lowest-energy 
minimum on the potential energy surface using global op- 
timization tools. The global minimum must be the equi- 
librium structure at zero temperature, and it is usually 
hoped that it will provide a good guide to the structure at 
higher temperatures. By locating the global minimum as 
a function of size a detailed picture of the cluster struc- 
ture and the non-monotonic variation of cluster proper- 
ties can be obtained. This can then allow a detailed-|Com- 
parison with experiments on size-selected clustersa and 
mass spectroscopic studies where the abunda^^es provide 
an indicator of particularly stable clusters. BSa 

However, global optimization becomes increasingly dif- 
ficult as the size increases, because, of the exponential 
increase in the number of minima.ErQ Furthermore, opti- 
mization of some sizes may be particularly difficult be- 
cause the pateatial energy surface has an unfavourable 
topography.Otffi These effects limit the size for which 
global optimization is feasible. For example, for the sim- 
ple Lennard-Jones (LJ) potential, after extensive study 
and computational effort, it is likely that the true global 
minimum for all. L J .. cl usters up to 150 atoms have 
now been found. B'E3llinl3 However, the last new putative 
global minimum in this size range was found as recently 
as 1999. t3 Plausible putative-global minimum have also 
been found up to Af=309oO For more complex poten- 
tials the present maximum feasible size is of course likely 
to be considerably smaller. 

Some clusters already exhibit bulk-like structures in 
the sizes accessible by global optipajzation (e.g. the al- 
kali halides with rocksalt structure&j) , however many do 



not. Extreme examples are provided by sodium clusters, 
which exhibit icosahedral structures up to at least 20 000 
atoms,E3 and by boron suboxide particles, which can itavc 
Mackay-like icosahedral structures with 10 14 atoms.li3 

Therefore, to go beyond the maximum size feasible 
for global optimization a different approach is necessary. 
This is possible if one settles for a more coarse-grained 
picture of the size evolution of cluster structure. So, 
rather than focussing on the the detailed size dependence, 
the general competition between structural types is in- 
stead followed. The question is then what type of struc- 
ture is most likely to be lowest in energy in a particular 
size regime, rather than what is the particular arrange- 
ment of the atoms in the lowest-energy minimum. 

Typically the energies of particularly stable sequences 
are compared to locate the crossover size at which the 
structure of the lowest-energy sequence changes. This 
method is commonly used an d s O |. cr ossover sizes have 
been estimated for LJ clusters JiSnJul metal clusters fer 
elements suctuas tungstenj23 nickel,E2l spdium,cil lead,c3, 
iron,c3 gold r 7 ! calcium and|-s.trontiumJ23 rhodium and 
palladium,Ej and aluminium EH and molecular clusters for 
molecules such as N 2 J^I CO^ and SF 6 B 

To get accurate crossover sizes it is important that 
the lowest-energy sequences for each structural type are 
used. The shapes of the optimal sequences are those 
that most closely resemble the polyhedron derived from 
the Wulff construction (or .a_modified form of it for non- 
crystallographic structuresEj) , which minimizes the sur- 
face free energy. However too many of the st 
this approach compare sub-optimal clusters.tr 
For example, it is particularly common to compare Ino 
decahedra and fee cuboctahedra to Mackay icosahedra, 
because they can have exactly the same sizes, even 
though the most stable sequences are likely to more 
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closely resemble truncated octahedra and MarksE-3 deca- 
hedra. 

As only the energies of the sequences are computed, 
the crossover sizes derived are only strictly valid for zero 
temperature. Usually the possible effects of temperature 
on these crossovers are neglected, and it is hoped that 
they will provide a reasonable guide for higher temper- 
atures. This neglect is probably because a simple way 
to compute the entropic contribution to the free energy 
of a structural type has not been available. Similarly in 
experiments on the evolution of cluster structure, size is 
usually the only variable that is considered, and so a sin- 
gle (presumaWj^temnerature independent) crossover size 
is presentedJSME3'E3 

However, there is an increasing body of theoretical ev- 
idence that thermal effects can play a key role in deter- 
mining a cluster's structure. There are a growing num- 
ber of examples of solid-solid transitions in clusters where 
the structure changes from fee or decahedral to icosahe- 
dral as the temperature increases. Such transitions have 
been found to occur for those small U^uateijs_that have 
a non-ic«pliedral global minimum J1J'^E3cjH for gold 
clustersP 2 ^ 3 !, for clusters interacting with a variable range 
Morse potentialJptia and even for models of clusters of 
Ceo molecules. e90 For molecular clusters, where low- 
pressure solid-solid or orientational phase transitions oc- 
cur for bulk, there has also been much interest in the .sizf?-. 
dependence of the transition temperature, e.g. SFgS'C^I 

Here we examine in detail the effects of entropy on the 
equilibrium structure of LJ clusters. In Section |lj we 
focus on the solid-solid transitions in those LJ clusters 
that have a non-icosahedral global minimum. In particu- 
lar we examine the different contributions to the entropy 
in detail. Using the principles learned from these exam- 
ples, in Section III we introduce methods to estimate the 
free energy of stable sequences of structures. Thus we 
are able to construct structural phase diagrams showing 
how the structure depends upon both size and temper- 
ature. There has been a previous attempt to calculate 
such structural phase diagrams based upon macroscopic 
estimates of the. different contributions to the free energy 
of a structured However, it is unclear how valid the ap- 
proximations made in that study arc, in particular the 
applicability of the macroscopic concepts to a small clus- 
ter. Clearly, a much more satisfactory approach is, as 
here, to compute the free energies directly from micro- 
scopic properties. Finally in Section [V we consider some 
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of the other factors that determine cluster structure and 
how these effect the comparison with experiment. A-brief 
report of some of this work has already appearedEj 
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FIG. 1: Structures of the global minimum (upper half) 
and the lowest-energy icosahedral minimum (lower half) 
for those LJ clusters with a non-icosahedral global min- 
imum. The labels denote the size and energetic rank of 
the minima. 



II. SOLID-SOLID TRANSITIONS 

The atpins in the clusters interact via a Lennard- Jones 
potentialEil and so the total potential energy is: 
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where e is the pair well depth and 2 1 / 6 <r is the equilibrium 
pair separation. 

Most LJ global minima for N^-< 309 have structures 
based upon the series of MackayEi icosahedraBEJ How- 
ever, there are a number of non-icosahedral global min- 
ima. These occur at sizes where the non-icosahedral mor- 
phology has a near optimum shape, whereas the icosa- 
hedral structures have an incomplete outer layer. For 
N < 147 there are 8 such cases, whiph.-CQixespond to 



the fee truncated octahedron at A 




the Marks 
and the Leary 



decahedra at iV=75-77_and 102-104,1! 
tetrahedron at iV=98.Ej These global minima and the 
lowest-energy icosahedral minima for these sizes are de- 
picted in Figure [j]. In the range 147 < N < 309 there are 
a further eight non-icosahedral global minima which are 
clustered in two sets aratLad the complete Marks decahe- 
dra at iV=192 and 238 Mia 

For these examples it has been found that the global 
minima are not the most stable structure up until the 
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cluster melts, but that a low-temperature solid-solid 
transition occurs at which icosahedr al .. s t ructur es be- 
come thermodynamically more stable.lliHl3c3c3'cll How- 
ever, these transitions are difficult to probe because there 
are lagSA (free) energy barriers between the structural 
types,0O making the dynamics of the transition ex- 
tremely slow. They occur on time scales way beyond 
those accessible with conventional molecular dynamics. 
For example, at the centre of the transition for LJ38 the 
rate constant has been estimated to be 43s -1 (using pa- 
rameters appropriate for Ar).E3 The equivalent rates for 
the larger clusters are likely to be orders of magnitude 
smaller. 

The only simulation method that has been able to di- 
rectly compute the equilibrium thermodynamics of the 
solid-solid-transitions of any p£ these clusters is parallel 
temperingEj (and its variantsB). Exchange of configura- 
tions between Monte Carlo runs at different temperatures 
allows ergodicity to be achieved at the low temperatures 
relevant to the solid-solid transitions. However, even for 
the smallest example, LJ 38 , this method is computation- 
ally demanding E3'E2l ■— . 

The superposition methodEl'E3E3E3 provides an alter- 
native approach and one which is particularly useful for 
our present purposes because it allows us to analyse the 
different contributions to the entropy. In this approach 
the partition function is written as a sum over all the 
minima on the potential energy surface, i.e. Z = Zi, 
where Zi is the partition function for the basin of attrac- 
tion surrounding minimum i. As long as all relevant min- 
ima are represented, ergodicity is ensured. The method 
is particularly simple to apply to solid-solid transitions as 
the number of minima that contribute to the transition 
is comparatively few. The method can also be applied 
to the melting of clusters, but in this case a statistical 
representation of the distribution of minima has to be 
used, because the number of minima associated with the 
liquid is, for alLbut the smallest clusters, too large to 
be enumerated.EJ Another particularly useful feature of 
the method is that, by restricting the sum to a certain 
subset of the minima, the thermodynamic properties of a 
particular region of configuration space can be obtained, 
i.e. Za = ^2i£A^ii where Za is the partition function 
for region A. 

To apply this method to the solid-solid transitions in 
LJ clusters, we first generated large samples of low-energy 
minima for each cluster using established techniques for 
searching potential energy surf aces. EtoE3 Some of these 
samples had been generated previously in work on the 
characterization of the potential energy surface tc. 
raphy of LJ clusters using disconnectivity graphs! 
The samples included thousands of minima, which is far 
more than would actually significantly contribute to the 
low-temperature thermodynamics. We then have to as- 
sume a form for the partition function of each minimum. 
Here, we use the harmonic approximation, which in the 
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FIG. 2: Low temperature heat capacity peaks associated 
with the solid-solid transitions in those LJ clusters with 
non-icosahedral global minima. The curves are labelled 
by the size of the cluster. 



classical limit gives 



; exp (-(3Ei) 



(2) 



where 0=l/kT, k=3N — 6 is the number of vibrational 
degrees of freedom, Ei is the potential energy of minimum 
i, Vi is the geometric mean vibrational frequency, and 
rii=2N\/hi is the number of permutational isomers of i, 
where hi is the order of the point group. 

This harmonic superposition method (HSM) has been 
shown to give a correct qualitative picture offGlaster ther- 
modynamics (even of the melting transitiorHS) , which 
becomes increasingly accurate at lower tempera 
Accurate anharmonic forms for Zi are available J. 65 "'" 
but the derivation of the associated parameters charac- 
terizing the degree of anharmonicity adds an extra degree 
of complexity, which is unnecessary here given the low 
temperatures at which the solid-solid transitions occur. 

For L J38 the fee to icosahedral transition has been pre- 
viously shown to occur at roughly two-thirds of the melt- 
ing tempera turn and to give rise to a feature in the heat 
capacityocJOO'Ell For the other clusters the transitions 
give rise to peaks in the heat capacity at low temperature, 
as shown in Figure ||. For reference, melting transitions 
occur .ai-xpughly 0.2-0.3e/c -1 for LJ clusters in this size 
range.OL°J Using the definition Za = Zb at the transi- 
tion temperature, T ss , where A and B represent the two 
competing structural types, we have calculated T ss values 
for all the clusters (Table |). 

The transitions must clearly be driven by the greater 
entropy of the icosahedral structures, but there are a 
number of contributions to this entropy. Firstly, there is 
the configurational entropy due to the number of minima 
associated with a structural type. This has a component 
due to symmetry. The number of permutational isomers 



TABLE I: Estimates of the solid-solid transition temper- 
ature T ss and associated properties for those LJ clusters 
with less than 150 atoms that have non-icosahedral global 
minima. V 'a and T>b are the geometric mean vibrational 
frequencies and Ha and hs are the orders of the point 
groups of the lowest-energy non-icosahedral and icosahe- 
dral minima, respectively, ua and ng have been evalu- 
ated using Equation (0) at T ss 



N 


HSM 


T ss /ek 
Einstein 


i 

Eq. (| AE/e 




n B /n A 


h A /h 


38 


0.121 


0.199 


0.316 


0.676 


1.0200 


31.7 


4.8 


75 


0.082 


0.234 


0.119 


1.210 


1.0475 


91.8 


20 


76 


0.046 


0.223 


0.053 


0.510 


1.0446 


4.32 


2 


77 


0.048 


0.199 


0.057 


0.565 


1.0451 


6.56 


4 


98 


0.004 


0.009 


0.006 


0.022 


1.0135 


12.0 


12 


102 


0.013 


0.096 


0.014 


0.086 


1.0201 


2.00 


2 


103 


0.016 


0.116 


0.018 


0.107 


1.0204 


2.00 


2 


104 


0.007 


0.069 


0.008 


0.048 


1.0212 


2.00 


2 



decreases as the order of the point group of a minimum 
increases. This term particularly favours the icosahedral 
structures at iV=38, 75 and 98, because of the high sym- 
metry (Oh, D 5 h and Tj) of the non-icosahedral global 
minima. The other contribution to the configurational 
entropy comes from the number of distinct low-energy 
minima of a particular type. This term, of course, de- 
pends on temperature. At zero temperature only the 
lowest-energy minimum contributes to the partition func- 
tion, Za, but as the temperature increases other low- 
energy minima begin to become populated. Furthermore, 
the contributions of other minima are greater the closer 
they are in energy to that of the lowest-energy minimum 
of that type. This term again favours the icosahcdra 
because there are many low-energy icosahedral minima 
with different arrangements of the atoms in the incom- 
plete outer layer E 

Secondly, there is the vibrational entropy. For the 
current examples the icosahedral structures have a 
smaller mean vibrational frequency (Table @), which 
again favours the icosahedra. However, unlike the previ- 
ous two contributions to the entropy, this term favours 
the icosahedra for a LJ cluster of any size, as illustrated 
by Figure where V (normalized to remove the non- 
specific size-dependence) is plotted for the LJjv global 
minima. Although there is a considerable variation in 
77 for the icosahedral structuresO (the complete Mackay 
icosahedra are the most rigid and those with an incom- 
plete outer layer can have values considerably smaller 
than Vya), only the vibrational frequencies for the non- 
icosahedral global minimum are greater than 77mi- 

We can quantitatively analyse these contributions to 
the entropy and their effect on T ss . We can neglect 
the differences in vibrational entropy by applying an 
Einstein-like approximation, where we assume that all 
the minima have the same mean vibrational frequency. 
The resulting values for T ss are also given in Table |. 




FIG. 3: Geometric mean vibrational frequency of the 
LJ global minima for N < 150. For those sizes with a 
non-icosahedral global minimum the value of V for the 
lowest-energy icosahedral minimum is represented by a 
cross. The zero, 77mi, is a fit to the vibrational frequencies 
of the complete Mackay icosahedra using Equation (||) . 



Although the transitions still occur, they do so at signifi- 
cantly higher temperature with the error increasing with 
size. For the largest cluster the estimate of T ss is roughly 
ten times too large. This increasing importance of the 
vibrational entropy is because v is raised to the power 
3iV — 6 in the expression for the partition function of a 
minimum in Equation |^. 

We can quantify the effects of the configurational en- 
tropy, if we rewrite Za as 

n A exp(-(3E A ) , , 

(Phv A ) K 

where tia, the effective number of minima (both per- 
mutational and geometric isomers) associated with A, is 
given by 

m exp(-j3(£j - E A )) 



n A (T) = J2- 



(vijv A Y 



(4) 



Using the properties of the lowest-energy minimum in 
region A to define Ea and 17 a, we give the relative values 
of the effective number of minima in each region at the 
solid-solid transition in Table |j} We also give the values 
of Iia/Iib, which correspond to the values of Ub/ua at 
zero temperature, where only differences in symmetry 
contribute. 

As expected the effective number of icosahedral min- 
ima is greater. However, for the four largest clusters 
tib/tia = hA/hs at T ss because, although there are many 
low-energy icosahedral minima, the transition occurs at 
too low a temperature for them to be thermally popu- 
lated, and thus to contribute to rig. 
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FIG. 4: Scatter plots showing how the geometric mean vibrational frequency of a minimum depends upon its energy 
and structure for (a) LJ38, (b) LJ75, U )nL|J?)|< and (d) LJio2- Isomers with different structures have been differentiated 
using orientational order parameters The unit of frequency is (e/mcr 2 ) 1 / 2 . 



From Equation |^ we can obtain an expression for T 5 

AE 



k (log (ub/tia) + Klog (va/vb)) 



(5) 



where AE=Eb — Ea- Of course, this equation cannot be 
straightforwardly applied to find T sa because ub/tia is a 
function of temperature. However, if we ignore the differ- 
ences in configurational entropy, this expression reduces 
to 



T — 



AE 



knlog {va/vb) 



(6) 



The values of T as that results from applying this expres- 
sion are given in Table |[ The values become noticeably 
more accurate as the size increases. This is because the 
vibrational term increasingly dominates over the config- 
urational term in the denominator of Equation (^|). 

From this analysis we have shown that differences 
in vibrational entropy are crucial to the occurrence of 



these solid-solid transitions. This conclusion is further 
strengthened by the absence of any solid-solid transi- 
tions when the LJ global minimum is icosahedral. These 
results, therefore, suggest that solid-solid transitions, 
rather than being unusual, should be expected for sys- 
tems where different structural types have systematic dif- 
ferences in v at sizes where the morphology with smaller 
vibrational entropy is the lowest in energy. 



III. STRUCTURAL PHASE DIAGRAMS 

Our aim is to produce a structural phase diagram 
showing how the LJ cluster structure depends on both 
size and temperature which is applicable to large clus- 
ters. At the relevant sizes (the energetic crossovers occur 
for clusters with thousands of atomsE^l) it is clearly not 
feasible to obtain large samples of minima as we have 
done for the examples considered in the previous section. 
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Methods are available to estimate.-the energy distribu- 
tions of minima for large clusters,c3 however these are 
involved and although possible for individual clusters are 
not feasible for a wide range of sizes. 

Equation (||) potentially provides a more readily ap- 
plicable approach, and as we have seen from the last 
section, this expression becomes more accurate as the 
size increases. Furthermore, its neglect of the differences 
in configurational entropy becomes more realistic in the 
coarse-grained approach we are seeking. Although there 
are likely to be significant differences in the numbers of 
low-energy minima for morphologies A and B for many 
sizes (as for some of the examples in the previous section) , 
for some sizes ha > n E but for other sizes n E > nA and 
so on average they will be of the same order. 

Equation also ignores the potential variation of v 
for minima of the same structural type. Therefore, if we 
are to use Equation (0), this variation should be signifi- 
cantly smaller than the difference in V between structural 
types. Figure [| illustrates that this is true using a selec- 
tion of the examples from the last section. The geometric 
mean frequency is a good order parameter to distinguish 
the structure of minima, and the two bands of minima 
become increasingly distinct as the size increases. 

It is also noticeable that for both structural types in 
each diagram there is a trend for the vibrational frequen- 
cies of the minima to linearly decrease as the energy 
increases. This is not unexpected. The higher energy 
minima tend to involve a less compact arrangement of 
the surface with fewer nearest-neighbour contacts and so 
have a lower frequency. The greater vibrational entropy 
of these minima is part of the reason why they are likely 
to become populated as the temperature increases. The 
creation-jof .surface defects and increased surface diffusion 
result EJoE3 However, as the trends in V with energy are 
similar for minima of both structural types, the use of 
Equation (|J) is not compromised. 

To use Equation (^|) to determine how crossover sizes 
depend on temperature we need expressions for the size- 
dependence of AE and F^/zTg. Of course, it would be 
unfeasible to follow the detailed non-monotonic size vari- 
ation of these properties over the size range we are inter- 
ested in, because this would require the global optimiza- 
tion of many very large clusters. Instead, we can use the 
properties of stable sequences of structures to construct 
a coarse-grained structural phase diagram. 

The energies and vibrational frequencies of these stable 
sequences can be fitted to the forms 

E{N)=a E N + b E N 2/3 + c E N 1/3 + d E (7) 

where the first two terms represent volume and surface 
contributions, respectively. These expressions can then 
be input into Eq. (|J) to map out the solid-solid transi- 
tions in the structural phase diagram. To complete the 
diagram, the melting line has to be determined. Again 
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FIG. 5: (a) Energies of the competing structural types 
for LJ clusters, (b) Mean vibrational frequencies for 
Mackay icosahedra, Marks decahedra with grooves of 
depth one and fee truncated octahedra with regular 
hexagonal {111} faces. In both parts, the data points 
represent specific clusters, the continuous lines are fits 
to the data using Eq. (^) and (||), and the zeroes, -Etoct 
and Ftoct, are a fit of the respective property for the fee 
truncated octahedra with regular hexagonal {111} faces. 



it will have the form 
T m (N) = a. 



dm 



(9) 



N 1 / 3 N 2 / 3 

The structural types that we should consider for LJ 
clusters are icosahedral, decahedral and fee. We have not 
considered unstrained close-packed clusters other than 
fee. Generally, fee clusters will have a lower surface en- 
ergy than hexagonal close-packed (hep) clusters because 
it is easier to form fee structures with a higher propor- 
tion of {111} faces. Therefore fee clusters are usually 
favoured over hep. However, a E is slightly lower for hep 
rather than fee clusters due to a more favourable next- 
nearest neighbour shell, and so for sufficiently large sizes 
hep clusters should be lower in energy. Similarly, it can 
be energetically favourable to introduce twin planes into 
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FIG. 6: Examples of large icosahedral, decahedral and 
fee clusters with optimal shape. The sizes are as labelled. 



fee clusters The differences in energy, though, are very 
small and this behaviour is different from the rare gases 
which form fee crystals. Therefore, we do not attempt to 
model this subtle structural competition. 

The energies of the most stable sequences for icosahe- 
dral, decahedral and fee clusters are represented in Fig- 
ure ^a, with the most stable sequences based on Mackay 
icosahedra, Marks decahedra and truncated octahedra, 
respectively. However, for each structural type there are 
small changes in the optimal shape as the size increases, 
as forms that more closely ijesemble the (modified) Wulff 
polyhedra become possible. El This gives rise to the lobes 
in each of the lines. By A=549 it becomes favourable 
for the vertices of the Mackay icosahedra to be unoccu- 
pied. As the size increases, {110} facets of increasing size 
are introduced to make the Marks decahedral and trun- 
cated octahedral forms more rounded,_and the grooves in 
the Marks decahedra become deeper .EiJ Three examples 
of structures that lie on the optimal lines are depicted in 
Figure |^ 

Up to TVwlO 000 the energies of these stable sequences 
of structures were obtained by minimization of all de- 
grees of freedom. Above this size minimizations were 
possible up to A«35 000 for decahedral and A«80 000 
for fee structures if the cluster was constrained to main- 
tain the correct point group symmetry [D^h f° r deca- 
hedral and Oh for fee). This constraint can reduce the 
computational cost by up to hi/2 (i.e. 60 for Ih, 24 for 
Oh and 10 for D§h) while the errors in the energies of the 
minima arising from the use of symmetry are negligible. 
These sizes are sigaificajitly larger than have been used 
in previous studies .E§E3'E!l 

As the energetic crossover between decahedral and fee 
structures lies beyond the largest size for which we could 
minimize the energy, we also had to use Equation ([?]) 
to obtain fits to the low-energy envelope formed by the 
succession of optimal sequences, which would be valid to 
this crossover region. The resultant energetic crossover 
sizes (Table || are consistent with or an improvement 
upon the most comprehensive of the previous studies.EH 

In order to calculate v the Hessian matrix must be di- 
agonalized. As the Hessian is not sparse — no cutoffs are 
applied to the potential so all atoms interact with each 
other — the cost of this operation scales as iV 3 . This scal- 
ing limits the sizes for which v can be calculated. It is un- 



TABLE II: The crossover sizes at zero temperature and 
at the melting point for quantum xenon, argon and neon 
clusters, modelled using an appropriate value of the de 
Boer parameter A, as compared to classical Lennard- 
Jones clusters. 







Aicos 


— »deca 


Adcca 


— >fcc 




A 


T = 


T = T m 


T = 


T = T m 


LJ 


0.000 


1690 


7440 


213 000 


6 670 000 


Xe 


0.010 


1810 


6510 


269 000 


5 520 000 


Ar 


0.030 


2130 


5500 


439 000 


4 380 000 


Ne 


0.095 


4640 


6230 


2 844 000 


6 390 000 



likely that continuum elastic theory can provide a short- 
cut to obtaining the required vibrational information for 
large clusters, since, as a recent study of two-dimensional 
clusters showed, it only begins to give a good description 
of the lowest-frequency modes when the size is roughly 
30 atomic diameters £3 

We have been able to characterize the vibrational prop- 
erties of the three sequences shown in Figure |5|b up to 
about 10 000 atoms, and other sequences up to about 
3500 atoms. As a v represents the volume contribution 
to the mean vibrational frequency, its value should be 
the same for all sequences of a particular structural type. 
Therefore, we used our most well-characterized sequences 
to determine a„ for each structural type. This gave 



3 =0.50970, 



l =0.51930, 



~fcc 



=0.52020. 
value, 



For other 
and then 



sequences we used the appropriate a L 
fitted the three remaining parameters. 

To construct the melting line a m is assigned the-jvalue 
of the zero pressure bulk melting temperatureLa and 
the other three parameters in Eq. (@) are fitted us- 
ing the melting points of the first four Mackay icosa- 
hedra, which were obtained from Monte Carlo simula- 
tions. For the 309-atom cluster parallel tempering proved 
to be necessary to ensure ergodicity. For this cluster a 
low-temperature shoulder was seen in the heat-capacity, 
which was indicative of a surface transition.l 76 H 7S l 

The structural phase diagram that results from these 
calculations is shown in Fig. |^. The phase bound- 
aries divide the plane into regions where the major- 
ity of clusters have a particular equilibrium structure. 
As the size increases, there is the expected progression 
from icosahedral to decahedral to fee clusters. The less- 
strained structures become favoured at larger N as the 
effect of their lower ab values dominates over their less 
favourable surface energies. For the sequences of Figure 



3 =-8.5320, 



=a =-8.6027, ag c =-8.6101. Fur- 
thermore, the effect of the vibrational entropy can be 
clearly seen from the slopes of the phase boundaries in 
Fig. ffl and from the differences between the crossover 
sizes at zero temperature and the melting point (Table 
||). At higher temperatures icosahedra and Marks deca- 
hedra remain most stable up to considerably larger sizes 
than would be expected from their energetic crossovers. 
This is because of the relative values of their vibrational 
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FIG. 7: Structural phase diagram for LJ clusters. The 
data points represent the melting temperatures of the 
four smallest Mackay icosahedra. 



frequencies: a^ os and aj; cca are 2.06% and 0.19% less than 
a l ^ c , respectively. Lower values of a„ are associated with 
structural types with a greater strain energy; the physical 
origin of this correlation will be explored elsewhere. 

Interestingly, the effect of entropy is greater for the 
decahedral to fee transition, even though the difference in 
v is smaller. The larger value of N in the denominator of 
Eq. (||) and the much smaller difference in more than 
compensate for the smaller difference in V. For example, 
it can be easily shown that 



dN 



(T = 0) 



Aa E + 2Ab E /3N^ 3 + Ac E /3N' 



T2/3 



kK\og(v A /v B ) 



Aa E 



kn\og{u A /v B ) 



for large TV, (10) 



where Aa E =a\ 

There has been much interest in the phase changes in 
clusters as the finite-size analogues of bulk phase transi- 
tions. The finite size, which can lead to unusual features, 
such as negative heat capacities, EElEll causes any transi- 
tion to occur not at a single temperature, but instead, 
both phases can coexist over a range of temperature.E3 
This is illustrated by the finite width of the peaks in Fig- 
ure ||. We can easily estimate this coexistence range for 
our examples if we define the upper and lower limit of 
this coexistence range, T^, as the temperatures at which 
(p A ,p B )=(0.1,0.9) and (0.9,0.1). This gives 



AE 



kn\og(v A /v B ) T log 9 



(11) 



However, as the size at which these transitions occur is 
large, the log 9 term does not significantly effect the de- 
nominator and so the coexistence range is very small. 
When T± are plotted on Figure the lines are indistin- 
guishable from those for the mid-point of the transition. 



The assumption of harmonicity was one of the approx- 
imations used to construct the structural phase diagram 
of Figure 0. It has been previously found that this ap- 
proximation leads to an overestimate of roughly 10% in 
the value of the melting temperature for small LJ clus- 
ters, because-,thje liquid-like state is significantly more 
anharmonicSEa One would expect the errors to be less 
for T ss because the differences in anharmonicity between 
two solid forms are likely to be much less and because 
the relevant temperatures are often much lower. 

One can introduce the effects of anharmonicity via 
temperature dependent frequencies, i.e. v(T) = i 
/3°//3), where j3° is a measure of the anharmonicity.l 
This then gives for the transition temperature: 



1 



K (log(^K)+log(^g 



IP 



kT: 



anharm 



AE 



(12) 



As the temperature occurs on both sides of this equation, 
it would have to be solved self-consistently. Alternatively 
an approximate solution can be obtained by taking the 
first terms of the binomial expansion of the second log 
and then substituting in the harmonic expression for T ss . 
This gives 



1 



Klog(P°/4) , P°b-P°a 



j^'J 1 anharm 



AE 



(13) 



The second term represents an anharmonic correction to 
^jharm _ rp Q a ppjy this equation, estimates of the size de- 
pendence of /3° would be required, which could then be 
fitted to a form equivalent to Equation (||) . Methods are 
available to obtain /3oJ£3 but unfortunately we were not 
able to extend these to sufficiently large sizes to enable 
fits of sufficient accuracy to be generated that could be 
applied to the large sizes relevant to Figure ^. 

We have also so far assumed that the clusters behave 
classically. However, it has recently been shown that 
quantum effects can have a significant effect on the ther- 
modynamics of LJ clusters.tiH We can apply these same 
methods to obtain quantum structural phase diagrams, 
which are applicable to the rare gases. The quantum 
equivalent of Equation ([|) is 



Za{T) 



n A exp( 



nr i - expR^r 



(14) 



where E A = E A + nh(v A )/2 includes the zero-point en- 
ergy, and (v A ) is the arithmetic mean frequency. To ap- 
ply this equation to find T ss it needs to be further sim- 
plified since it would require the characterization of the 
size evolution of the whole frequency distribution. To il- 
lustrate the potential effects of quantum derealization, 
we choose the simplest approximation, namely to repre- 
sent the distribution by n modes at the mean frequency. 
Ignoring differences in the configurational entropy as be- 
fore, this then gives 



AE° 



knlog (j^ 



cxp(— [jhV B ) 



(15) 
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FIG. 8: Structural phase diagrams for rare gas clusters, 
illustrating the effect of the introduction of quantum ef- 
fects. The classical results are solid lines (LJ), and the 
phase boundaries for the rare gas clusters are as labelled. 



where AE° = E% - E° A . 

The introduction of quantum effects can also lead 
to significant changes in the melting temperature of 
clusters. Ellc^l Following Ref. ^ we simply model this 
change by 

T%(N,A) = T rn (N) - A 2 (a Q + -^j (16) 

where the coefficients oq = 5.18 and 6q =j_6.84 have been 
obtained using results for small clusters,E3 and A is the 
de Boer parameter. A is a measure of the quantum de- 
localization and corresponds to the value of h in reduced 
units, i.e. A = Ti/a^/me. 

We have calculated structural phase diagrams for val- 
ues of the de Boer parameter that are appropriate to 
some of the rare gases. The results are illustrated in 
Figure || and Table |l[ There are two competing effects 
on the structural phase diagrams resulting from the in- 
troduction of quantum delocalization. Firstly, the zero- 
point energy stabilizes those structures with lower fre- 
quencies, and therefore the zero temperature crossover 
sizes increase substantially as A increases (Table ||. In- 
deed the differences between neon and the other rare 
gases at low temperature are large enough to be poten- 
tially observable, if equilibrium clusters can be grown 
in this temperature range. Secondly, for temperatures 
lower than the vibrational temperature T V ib = hV/k, the 
vibrational contribution to the entropy is much less in 
the quantum case. This is a corollary of the fact that 
for a quantum harmonic well C v j (N — 2) is close to zero 
at low temperature and only rises near to 3 above T v ib, 
whereas classically C V /(N — 2) = 3 for all T. There- 
fore, the temperature dependence of the crossover sizes 
are diminished, and the phase boundaries initially rises 
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FIG. 9: A comparison of the energy interpolated be- 
tween particularly stable sequences and the energies of 
the global minima. 



vertically from the x-axis and in some cases even cross 
the classical line. 



IV. DISCUSSION 

The development of the structural phase diagrams pre- 
sented here provides the starting point for a more sophis- 
ticated discussion of the size evolution of cluster struc- 
ture, which emphasizes the effect of both size and tem- 
perature. The one-dimensional picture which focusses on 
the effect of size and neglects thermal effects should be 
discarded. Therefore, we hope that this work will pro- 
vide a stimulus for new experiments and calculations that 
seek to illuminate the role of temperature in determin- 
ing structure. Jndeed, recent experiments, on clusters of 
Ceo molcculesE-3 and strontium clustcrsE3 illustrate the 
interesting results that can be achieved. The methods 
developed here should be applicable to a wide range of 
materials, for which empirical potentials are available. 
An illustrative example fac silver clusters modelled by, 
the Sutton-Chen potential^ has already been provided.^ 
This system has a structural phase diagram which was 
similar in form to that for LJ clusters. However, it would 
be particularly interesting to map out the different pos- 
sible families of phase diagrams. For example, in sodium 
clusters there would also be additional competition be- 
tween the icosahedral, decahedral and fee structures con- 
sidered here and clusters that are body-centred cubic. 
These methods could also be potentially applied to clus- 
ters in a matrix or on a surface. 

This point of view should also help in comparisons 
of experiment and theory, and of different experiments. 
Firstly, it helps one not to jump to hasty conclusions. For 
example, if two experiments come to different structural 
conclusions for clusters of the same size, this does not 
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necessarily imply that there has been a misinterpretation 
of the data; the differences might simply reflect different 
temperatures. Similarly, the failure of energetic calcu- 
lations to reproduce the experimental observations does 
not necessarily imply the inadequacy of the theoretical 
description of the interatomic forces, because free ener- 
gies need to be calculated for a rigorous comparison. Sec- 
ondly, and more positively, theoretical and experimental 
agreement on the size and temperature dependence of 
structure would provide a strong vindication of the ap- 
proaches used. 

There are a number of aspects that should be born in 
mind when interpreting the structural phase diagrams. 
Firstly, they only provide a coarse-grained picture of the 
size evolution, dividing the phase diagram into regions 
where the majority of sizes have the same equilibrium 
structure. To construct the diagrams we have assumed a 
monotonic dependence of properties on N (Equations |t], 
|§| and ^). However, as is well-known, most cluster prop- 
erties show a complex non-monatonic size dependence, 
particularly for smaller clustersEH Indeed, it is these size 
effects that are often of most interest. 

This non-monotonicity is illustrated for V in Figure || 
and for E in Figure |^. Similarly, the variation m-T—. 
for small LJ clusters has been catalogued by Frantz.Ej'E2l 
It is noticeable that the variations in E and V can be 
of a similar magnitude to the systematic differences be- 
tween the structural types. For example, at three points 
the energies of the global minima in Figure [| coincide 
with those for truncated octahedra (N=38) and complete 
Marks decahedra (iV=75,192). Consequently, the real 
phase boundaries when resolved at the level of individ- 
ual sizes will be much rougher than in Figure [?]. Indeed, 
the examples in Section || correspond to small islands of 
stability for fee and decahedral structures, within a sea 
of icosahedral dominance at small N. 

The relative importance of these variations does dimin- 
ish with increasing size as the proportion of atoms in the 
surface decreases. Consequently, the coarse-grained ap- 
proach used for the structural phase diagrams becomes 
more reliable at larger sizes, but should be applied with 
caution when crossovers occur for sizes with less than a 
few hundred atoms. There are other implications of this 
decrease in fluctuations. For example, it is clear from 
Figure ^| that the 38-atom truncated octahedron will be 
the only fee global minimum up to at least iV=1000. 

Secondly, the position of the phase boundaries in the 
diagrams sensitively depends on the, usually small, dif- 
ferences in the energetic and vibrational properties, in 
particular Aa and Ab, between structural types. It is 
asking a lot of an empirical potential to reproduce these 
properties accurately enough to give phase diagrams that 
could agree in detail with experiments. For the current 
example, where different structural types are all based 
in some sense upon close-packing but perhaps involving 
multiple twinning and hence strain, the crucial ingredi- 
ents for the energetics would be a correct reproduction of 
the relative surface energies for {111}, {100} and {110} 



faces, and the elastic moduli that determine the energetic 
penalty associated with the strain. Similarly, one would 
also want the potential to reproduce not only the bulk 
phonon dispersion curves, but also the effects of surfaces 
and strains on the vibrational properties. Furthermore, if 
one were to apply this approach to clusters, where, as for 
sodium, two different bulk structures could be adopted, 
the potential would also have to be able to reproduce 
the properties of both forms. However, most empirical 
potentials are only fitted to a small subset of these prop- 
erties, usually those associated with the bulk. Therefore, 
given the current limitations associated with empirical 
potentials and the current impossibility of performing ab 
initio calculations for the relevant size ranges, it is per- 
haps better to use the tools developed here to understand 
the qualitative form of the structural phase diagram one 
would expect for a particular material. 

Thirdly, the experiments one is seeking to compare 
with may not be at equilibrium. Evidence for the non- 
equilibrium nature of many clusters has come from a 
number of sources. It is suggested by the mixture of 
structures often observed for metal clusters, for exam- 
ple, when .supported clusters are-probed by electron 
microscopyEa or x-ray diffractionpJ greater structural 
uniformity would be expected for samples of equilib- 
rium clusters. More rigorous evidence for the kinetic ori- 
gin of this diversity comes from experiments that have 
shown that the structural composition can bei-changed 
simply through the cluster source conditions pi'tB and 
that well-defined structures are only obtained after care- 
ful annealing.Ea Similarly, it has recently been shown 
for clusters of Cgo molecules that the equilibrium non- 
icosahedral structuies_are only located after annealing at 
high temperature 

This difficulty of achieving equilibrium is also sug- 
gested by a range of recent theoretical results. For exam- 
ple, the potential energy surfaces of the examples in Sec- 
tion [n| have a multiple-funnel character, where the differ^, 
ent funnels correspond to the different structural types.Q 
This multiple-funnel character is not restricted to these 
examples but appears to be ubiquitous .Ol3oE3 As the 
(free) energy barriers between funnels are typically large, 
the dynamics of structural transformations can be very 
slow, giving rise to the possibility of kinetic trapping in 
out-of-equilibrium structures. . .. . 

This effect is seen in recent growth simulations. 
Typically, above a certain size, which is dependent on 
the material, temperature and growth rate, structural 
transformations are unlikely, and so new growth then oc- 
curs around this preserved seed cluster. The results are 
particularly dramatic for (C6o)jv where perfect Mackay 
icosahedra are invariably formed, even though they are 
far from equilibrium.^ Similarly, it has been found that 
freezing oi, molten clusters give rise to a a mixture of 
products,B3 where the relative amounts reflect the width 
of the respective funnels JEj rather than the thermody- 
namics. These results emphasise that it will be diffi- 
cult to obtain equilibrium when a structural type is only 
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favoured at low temperature. This is particularly the case 
for the larger sizes often associated with the decahedral 
to fee crossover. 

Many of the above considerations are relevant when we 
seek to compare our results for Lennard- Jones clusters 
with experimental results for rare gas clusters. The first 
difficulty is that a clear-cut picture has not yet emerged 
from the experiments. The most easy to interpret data 
comes from mass spectroscopic studies of rare gas clus- 
ter ions. Clear magic numbers that are signatures of 
icosahedtai-packing have been obs erved f or small clusters 
of argon pi!23 krypton and xenon.l 101 H 102 l The icosahedra 
have been observed up to iV=923, with-^ao sign of any 
crossover to a different structural type.t 02 ] However, the 
situation is less clear for neon. In the only mass spec- 
troscopy experiment on neon clusters, of which w£— are 



aware, no strong signal of shell effects was found.liril It 
should also be remembered that there may be differences 
between the structure of ions and-jieutrals, as illustrated 
by recent calculations for Ne^,,ll2a although these differ- 
ences would be expected to be less important as the size 
is increased. 

There have also been a series of experiments by 
slectron diffraction to probe argon 
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To reach a structural assignment the 
experimental diffraction patterns have to be compared 
to those calculated for structures of the appropriate size. 
Although this provides good evidence that small clusters 
are icosahedral with a suggested crossover from icosahe- 
dra at JV ~ 750, it has proved to be much harder to 
obtain a good fit to experiment for larger clusters with 
neither decahedral or fee structures providing adequate 
agreement. Van de Waal has investigated a series of mod- 
els that provide better fits, .which include features, such 



as crossed stacking faults 
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crossed twins 
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that 



could lead to faster growth. For very large clusters, how- 
ever, adequate agreement could only be obtained from 
polycrystalline clusters-involving fee, hep and random 
close-packed domains.liiil These results again seem to 
confirm the difficulty of achieving equilibrium, and are 
further reinforced by simulations of the freezing of LJ 
clusters, in which a mixture of structural forms is ob- 



served for larger clusters. til 

Somewhat different conclusions have been reached 
from EXAFS experiments on argon clusters,E£l in which 
fee structures were suggested for clusters with more than 
200 atoms. However, this conclusion was based on com- 
parison with results obtained for a very limited set of 
candidate structures, and no estimate of the experimen- 
tal temperature was given. 

Given the above experimental situation it is not sur- 
prising that there is not yet full agreement between ex- 
periment and theory. Furthermore, although the LJ po- 
tential provides a reasonable description of the proper- 
ties 
curve 



of the structural phase diagram to be robust. There- 
fore, it would first be interesting to see if an icosahedral 
to decahedral crossover could be revealed for argon clus- 
ters as the experimental clusters are brought closer to 
equilibrium. Then, if this could be achieved, it would 
be interesting to see if the temperature dependence of 
the crossover size that we have predicted here was repro- 
duced. 

The results we have obtained for LJ clusters have high- 
lighted the important role that differences in vibrational 
entropy can play in determining structure. Striking ex- 
amples for the considerable effect of vibrational proper- 
ties on phase diagrams have also been q btain 
thermodynamics of substitutional alloys 
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larly, in the field of supercooled liquids these 

tionEli; 



pr the 
B3 Simi- 
have been 

recent suggestions, based upon simulation! 11 ^ and calcu- 
lations for model potential energy surfaces ] 119 1 that vi- 
brational properties could play a key role in the complex 
interplay that gives rise to the distinction between fragile 
and strong liquids. 
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